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Abstract We study the effects of synaptic plasticity on the determination of firing 
period and relative phases in a network of two oscillatory neurons coupled with re- 
ciprocal inhibition. We combine the phase response curves of the neurons with the 
short-term synaptic plasticity properties of the synapses to define Poincare maps for 
the activity of an oscillatory network. Fixed points of these maps correspond to the 
phase-locked modes of the network. These maps allow us to analyze the dependence 
of the resulting network activity on the properties of network components. Using a 
combination of analysis and simulations, we show how various parameters of the 
model affect the existence and stability of phase-locked solutions. We find conditions 
on the synaptic plasticity profiles and the phase response curves of the neurons for 
the network to be able to maintain a constant firing period, while varying the phase 
of locking between the neurons or vice versa. A generalization to cobwebbing for 
two-dimensional maps is also discussed. 
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1 Introduction 

The output of a neuronal network, determined in part by the relative spiking times of 
its individual neurons, depends on the coordinated activity of its neurons. Observed 
phase relationships result from the combined effects of individual cells and synaptic 
connections whose properties change dynamically. For example, individual neurons 
in a network can differ in their intrinsic properties, being silent, spiking or bursting; 
different neurons can have different responses to the synaptic inputs they receive, 
and the synaptic inputs themselves can differ widely. These different characteristics 
all play a role in determining the resulting network activity. Determining how these 
dynamically varying components work together to influence the network activity is a 
question of considerable interest. 

Many studies have explored the question of how period and phase are determined 
in an oscillatory neuronal network [1-7]. One of the main tools used in these studies 
is the phase response (or resetting) curve (PRC) of an individual neuron. The PRC 
measures how the phase of firing of an oscillatory neuron changes as a function of 
perturbations that it receives at different phases of its oscillation. In prior work [8], 
the PRC has been used to define a ID map that measures the degree of network 
synchrony. This map allows for the analysis of the network activity in a reduced 
system by considering only the effect of the synaptic inputs on cycle length, rather 
than considering multiple dynamic variables. Several studies used similar methods to 
study the activity of neuronal networks [2, 8-12]. PRC -based maps were also used 
to incorporate some properties of neurons or synapses. This approach was applied 
to understand synchronization of adapting neurons [2, 5] as well as the effect of 
conduction delays on network synchrony [1, 13]. 

In the current study, we are interested in predicting phase-locking by deriving 
maps that combine PRCs with information arising directly from synapses that dis- 
play frequency-dependent short-term plasticity. Increases in presynaptic firing fre- 
quency can strengthen (facilitation) or weaken (depression) a synapse [14]. Some 
synapses show a combination of both, in which case the maximal synaptic amplitude 
is achieved at a specific presynaptic frequency [15] referred to as the preferred fre- 
quency of the synapse. Synaptic plasticity can be described with models having two 
variables, one for depression and the other for facilitation [15, 16]. 

The main advance in our work is the derivation of tools for analyzing higher- 
dimensional maps that incorporate the effects of synaptic plasticity and provide pre- 
dictions on circumstances under which an oscillatory network of neurons will phase- 
lock and at what period. In particular, we consider a network of two neurons, mutually 
coupled by inhibition in which the synaptic strength is frequency dependent. In deriv- 
ing these maps, we must not only track the phases of each cell, but also the strength 
of each synapse. As a result, the ID map that sufficed in prior studies needs to be 
replaced with 2D or 3D maps. For 2D maps, we derive a geometric method that gen- 
eralizes the idea of cob webbing. Namely, we show how iterations of the map can be 
tracked through different 2D surfaces. Moreover, projections of these surfaces onto a 
common plane yields two curves whose intersection is a fixed point of the map that 
corresponds to a phase-locked solution. We derive conditions on the PRCs and the 
parameters that govern synaptic plasticity of the neurons to show how a network can 
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have a range of parameters over which the network period remains constant, but the 
phase of locking between cells changes, or vice versa. We also show that the methods 
derived apply to networks that are heterogeneous either in the intrinsic properties of 
individual cells, in their synapses, or both. 

2 Model and Methods 

2. 1 Dynamics of Neurons 

We use Morris-Lecar (M-L) model neurons to conduct our analysis [17]. An iso- 
lated M-L neuron is modeled by leak (L), potassium (K) and calcium (Ca) current. 
The K current is driven by a dynamic activation variable w, while the Ca current 
depends on an instantaneous function ffi oo of the membrane voltage (V). When the 
neuron is synaptically coupled to another neuron, the synaptic current received from 
the other cell is also included in the equation governing the membrane voltage. For 
two M-L neurons coupled with synaptic inhibition, the equations for voltage V and 
K activation variable w are given by 



C—^ = /app - {gL(Vi - El) + gKwiVi - Ek) + ~gCamM){yi - Eca) 



for ij = 1,2, / ^ J, where m^(V) = 0.5(1 + tanh((y - VJ/V^)), w^(V) = 
0.5(1 +tanh((y - Vc)/Vd)) and r^(y) = l/(0cosh((y - Vc)/2Vd)). The conduc- 
tances (in nS) are gi = 2, gx = ^, gca = 4, the reversal potentials (in mV) are 
El = —60, Ek = —84, Eca = 120 for the leak, potassium and calcium currents, 
respectively. The synaptic reversal potential £'syn is —80 mV, modeling an inhibitory 
synapse. Due to the presence of the Heaviside function H(y — Vth), the synapses are 
all-or-none and activate (deactivate) instantaneously when the presynaptic voltage is 
above (below) the synaptic threshold Vth = 0 mV. The results described in this study 
remain qualitatively similar if the value of Vth is changed. 

Below we will provide more details of the synaptic conductance ^pre^post- We 
change the applied current /app (in pA) between 41.2 and 44.9 to obtain a set of 
intrinsic periods (in ms) ranging between 100.3 and 180.83. The rest of the model 
parameters are C = 20 pF, 0 = 0.067 (dimensionless), and Va = — 1.2, = 18, 
Vc = 12, = 17.4 in mV. Throughout the paper, units for time are in msec. 

2.2 Phase Response Curves 

The phase response curve (PRC) of an oscillator describes how the time of the next 
spike of an oscillator changes depending on the phase at which it receives a per- 
turbation (Fig. la). In general, the PRC can be computed numerically (for model 
neurons) or experimentally (for biological neurons) by injecting a brief perturbing 



Woo(Vi) - Wj 



Vth) ■ (Vi - £syn)) 



(2.1) 



dwi 
dt 
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Fig. 1 PRC due to synaptic input, a A brief perturbing current pulse stimulus (arrow) is used to measure 
the PRC as described in Eq. (2.2). b The PRC obtained from the Morris-Lecar model (2.1) neurons by 
inhibitory synaptic input. The parameters are /app = 42.2, synaptic conductance ^pre^post = 0.1 and 
synaptic duration ta = 14.3 

current (such as a small current pulse) and measuring the effect of this perturbation 
on the cycle length as a function of the phase of the perturbing input. If the pertur- 
bation is infinitesimally small, then an infinitesimal phase response curve (iPRC) of 
the model neuron can be obtained by linearizing the governing differential equations 
about the limit cycle and solving the adjoint equation. In this work, we will use the 
term PRC to refer to responses calculated by direct perturbations, for example ones 
that imitate synaptic inputs. 

Denote by Pq the intrinsic period of a cell. Suppose a perturbation is given at 
time dt after the firing of the cell. This yields a phase (j) = dt / Pq of the perturbation. 
Denote by P the time between when a cell fires prior to a perturbation and the sub- 
sequent firing of the cell when a perturbation is given at phase 0. Then we define the 
PRC as 

Po- P 

Z(0) = (2.2) 

We have chosen parameters so that in the M-L model oscillations arise through 
a saddle node on invariant circle (SNIC) bifurcation. Neurons that oscillate through 
a SNIC bifurcation have a Type 1 iPRC [18], which is always of one sign. In the 
case of an inhibitory perturbation received by the neuron, the Type 1 iPRC is never 
positive and the next firing time is therefore delayed. A PRC obtained from our model 
neurons for a specific synaptic strength is shown in Fig. 1. It is computed by applying 
a perturbation of the form 

^syn = ^pre^postH(Vpre ~ ^th)(^post ~ -^^syn) 

The reference point to compute the PRC is chosen to be when V crosses Vth in 
the positive direction. Note again that this method of computing the PRC is different 
from computing the iPRC of a spiking neuron which yields a strictly Type 1 PRC. 
The PRC we obtain is very similar, but there is a region of the PRC that is positive 
near small stimulus phases due to the longer active duration of the neuron. 
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2.2.7 Selection of PRCs 

In order for our analytical estimates to match the results of numerical simulations of 
the model, we took advantage of the computability of a PRC for the M-L neuron. In 
each iteration, we numerically computed the response of a neuron to a synaptic input 
of a specific strength at a specific phase. Although this method yields accurate results, 
it is computationally slow and it is almost impossible to implement on biological 
neurons. For this purpose, we created a meshed PRC measured at discrete phase 
points and for a discrete set of predetermined synaptic strengths. We used mesh sizes 
of 0.1 for the phase and 0.0125 for the synaptic strength to obtain a total of 77 points 
of numerically computed phase response values. The responses to the phases and 
strengths not on the mesh points were calculated by linear interpolation. 

2.3 Model for Synaptic Plasticity 

The short-term synaptic plasticity in spiking cells can be described by a phenomeno- 
logical model [16]. We modify this model for neurons that have broader action poten- 
tials or those for which the burst envelope instead of individual spikes are modeled. 
To account for the longer durations that the neurons spend above the threshold we as- 
sume that there are two variables which determine the strength of the synapses when a 
neuron fires; the depression variable (r) and a facilitation variable {u). The depression 
variable r represents the amount of available synaptic resources, while the variable 
u represents the amount of utilized synaptic resources. They change according to the 
activity of the presynaptic cell and together determine the synaptic strength. These 
variables obey the following dynamics: 

V< Vth 

(2.3) 

v>yth 
v<yth 

When the membrane voltage of the presynaptic cell is above the synaptic thresh- 
old Vth, the depression variable r approaches 0 with time constant ri, representing 
the depletion of available synaptic resources. During this time interval, the facilita- 
tion variable u approaches 1 with time constant rs representing the increase in uti- 
lized resources. When the membrane voltage is below the synaptic threshold, these 
variables recover to their steady-state values of 1 and U , with time constants T2 and 
r4, respectively. The strength of the synapses is determined by scaling the maxi- 
mal synaptic conductance by the product of the values of these variables when the 
presynaptic cell crosses Vth- If the presynaptic cell fires a sequence of spikes, then 
the term ^th cycle refers to the time duration between the ^th and n + 1st cross- 
ings of Vth- Hence, the synaptic conductance at the start of the n\h cycle is given by 
^pre^post = ^pre^post^w^w. where and Un are the values of r and u when the presy- 
naptic membrane potential passes synaptic threshold in the n\h cycle {n is defined 
below). 
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2.3.1 Steady State Synaptic Plasticity Profiles 

If the presynaptic cell fires with a fixed frequency, then it reaches an oscillatory steady 
state. The values r and u then also reach steady states and each oscillates between a 
minimum and a maximum value. At steady state, when crossing the synaptic thresh- 
old, Vn is at a maximum, rmax, while Un attains its minimum, i/min- These values can 
be calculated from (2.3) as 

1 _ e-^b/-c2 
^"^^^ " l-e-ta/rie-tb/n 

(2.4) 

^min — \ f /_ TTTV. 

where ta and t\) are the durations that the cell spends above and below Vth, respec- 
tively. 

It is often possible to measure the strength of the synaptic output when the presy- 
naptic neuron is driven in a range of frequencies. The values rmax and i/min depend on 
the presynaptic frequency and an appropriate choice of time constants allows for our 
model to fit a variety of frequency-dependent synaptic outputs. In particular, we are 
interested in synapses whose strength is maximal at a unique "preferred" frequency 
as we have observed in experimental measurements [19]. In our results presented 
below, we will use period instead of frequency for ease of analysis. By choosing ap- 
propriate parameters, therefore, we can match the period at which the peak of the 
product rmaxWmin is maximized with the experimentally measured preferred period of 
the synapse. We define the function 

^(/^)=^rniax(/^)^min(/^) (2.5) 

as the synaptic strength at the time of firing of a presynaptic neuron with constant 
period P = ta+tt,. We will assume that the changes in period of the bursting neurons 
affect only the inter-burst duration (i.e., ta is fixed). We will henceforth refer to this 
relationship (2.5) as the steady- state synaptic plasticity profile. 

Fi gure 2 shows plots of the steady-state values of rmax? ^min and the full synaptic 
plasticity profile (rmaxWmin) of a synapse as a function of the firing period, for a given 
set of parameters. Here = 15. The peak of the synaptic plasticity profile in this case 
occurs at P = 170. For ease of analysis in Sect. 3.3 and thereafter, we use a Gaussian 
function approximation for the steady-state synaptic plasticity profile giP): 

g(P) = 0.75^-^^-^p-f^'/^^^'^ + 0.75 (2.6) 

where Ppref is the peak of the profile corresponding to the preferred period of the 
synapse and a determines the spread. 

3 Results 

We derive Poincare maps that relate the firing times of a network of two neurons cou- 
pled with reciprocal inhibition. We assume a predetermined one-to-one firing order 
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Fig. 2 Steady-state values of 
plasticity variables. The 
maximum value rmax that the 
depression variable r and the 
minimum value u^i^i that the 
facilitation variable u reach at 
the steady state at the onset of 
presynaptic activity plotted 
against the presynaptic period. 
The plasticity profile of the 
synapse is given by their product 




^0 100 200 300 400 500 600 

between the neurons. The fixed points of these maps correspond to one-to-one firings 
of the neurons at the steady state. It is possible to derive similar maps assuming or- 
ders of firing that are not one-to-one, but these derivations are beyond the scope of 
the current study. We first assume a fixed synaptic strength between the neurons in 
Sect. 3.1. When the synapses have a fixed strength, only the phase response infor- 
mation of the neurons is used to determine the network activity, as has been shown 
previously [8]. In Sect. 3.2 we derive maps that describe the network activity when 
the synapses between the neurons are plastic. We compare two cases. In one case, 
we assume that the synapses obey the plasticity dynamics given in Eq. (2.3). In the 
second case, we consider synapses that obey the corresponding steady-state values 
given in Eq. (2.4). The latter case results in a lower-dimensional map. In Sect. 3.3, 
assuming both synapses obey a steady-state plasticity profile (2.6), we examine how 
changes in these profiles determine the network period and relative phase relations. 
We find conditions for a network to be able to keep a fixed firing period but vary the 
relative firing phase between its neurons, and vice versa. 

3.1 Map for Phase with Static Synaptic Strength 

We start with a network of two oscillatory neurons reciprocally inhibiting each other 
with constant synaptic strength. We will derive a ID map that measures the phase 
difference between the onset of firing of the two cells. A fixed point of the map 
corresponds to a 1 : 1 phase-locked solution. We then derive the criteria for existence 
and stability of fixed points. Finally, we test the map in a network of two M-L model 
neurons. 

Consider a network of two oscillatory cells, A and B, coupled with reciprocal 
inhibition (Fig. 3a). Assume that the synaptic strengths between the cells are constant 
in each spike, i.e., ^a^b = ^b^a = g- The intrinsic period of cell A and cell B are 
denoted by Pq and Qo, respectively. When the neurons are synaptically coupled, the 
time between subsequent firing of the same neuron may change. This time is called 
the cycle length, denoted by Pn and Qn in cycle n, respectively for A and B. 

We derive a Poincare map for the relative firing times of the neurons when they 
are synaptically connected. We choose the Poincare section to be at Va = Vth- The 



^ springer 



Page 8 of 29 



Z. Akcay et al. 



Fig. 3 Schematic diagram of 
the coupled network and the 
map variables, a Schematic of 
the network connectivity 
diagram, b The cycle length Pn 
of cell A in cycle n (measured 
for the M-L simulations when 
voltage crosses Vth) can be 
divided into the delay between 
cell A activity to cell B activity 
{dtn) and the opposite {dxn). 
The cycle period Qn of cell B in 
cycle n is dxn + dt^^i 





amount of time that passes after cell A fires until cell B fires is denoted by dtn, 
while the amount of time after cell B until cell A fires is denoted by dtn (Fig. 3b). 
The (activity) phase of neuron A (or B) is defined as the firing time dtn (or dtn) 
normalized by the cycle length. Therefore, the phases of A and B are, respectively, 
given by 

^n=dtn/Pn (3.1a) 
On=dTn/Qn (3.1b) 

In the derivations of the maps, we will make use of the PRCs of A and B which 
are defined in terms of Pq and Qo, the intrinsic periods of A and B. To simplify these 
derivations we introduce the notation of the "intrinsic phase" of neurons A and B 
which are defined, respectively, as 

(l)n=dtn/Po (3.2a) 

0n=dTn/Q0 (3.2b) 

We denote the PRC of cell A and cell B as Za(-) and Zb(-). respectively, for 
synaptic inputs with a fixed strength. Rewriting the PRC relationship (2.2) for the 
cycle lengths, we can obtain the cycle lengths of each cell in cycle n as 

Pn = Po{l-ZA(M) (3.3a) 

Qn = Qo{l-ZB(On)) (3.3b) 

The following equations relate the firing times of the two cells: 

dtn+dTn = Pn (3.4a) 

dxn+dtn+i = Qn (3.4b) 
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From (3.3a) and (3.4a), On can be written in terms of 0^: 



On = ^ = -^(Pn-dtn) = -^[Po{l " ZA(0n)) " /^O0n] 

Qo Qo Qo 

= — (1 - Za(0^) -0nj 
Qo 



(3.5) 



Similarly, 0^+i can be expressed in terms of On'. 

0/1+1 



dtn+l \ . ^ , . 1 



Po 

Go, 



(Gn - dxn) = ;^[Go(l - ZB(On)) " Go] 



= ^(1-Zb(^;,)-^^) 

Po 



(3.6) 



using (3.3b) and (3.4b). 

Thus, plugging Eq. (3.5) into Eq. (3.6) defines the following ID map for the intrin- 
sic phase of cell A (3.2a) when the 1 : 1 firing order between the cells is maintained: 



0^+1 = 77(0^) 

^ Go 
Po 



1-Zb(^^(1-Za(0J-0^)) -l + ZA(0n) + 0n (3.7) 



The condition for a 1 : 1 phase-locking solution is 0^ = 0^+i = 0*. Plugging this 
into the map gives the condition for a fixed point as 



Po(1-Za(0*)) = Go(1-Zb(^*)) 



(3.8) 



where 6** = ^(1 - Za(0*) - 0*). The fixed point is stable if |/7'(0*)| < 1, hence 
the stability condition is 



|(z^0*) + i)(z^(r) + i)|<i 



(3.9) 



This result was previously found in [8]. If the neurons are identical, Pq = Go and 
Za(-) = Zb(-) = Z(-). Then the map (3.7) reduces to 



0^+1 = 77(0^) 

= -Z(l - Z(0^) - 0^) + Z(0J + 0^ 

The fixed point equation (3.8) becomes 

Z(0*) = Z(1 -Z(0*) -0*) 

and the stability condition (3.9) becomes 

l(z'(0*) + i)(z'(i-z(r)-r) + i)|<i 



(3.10) 



(3.11) 
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In this symmetric case, the phase locking of the network does not depend on the 
intrinsic periods Pq of the network neurons. The phase of cell A (3.1a) in cycle n can 
be obtained from the relation 

~ dtn (pn Po 

which can be simplified using Eq. (3.3a) to 

0. = ^^^^/(0.) (3.12) 

1 - Z(0^) 

Given the map (3.10) for 0^, in order to derive a map for 0^+i, we need the 
function given in (3.12) to be invertible. The function / is monotone increasing in 
[0, 1] if and only if /^(0) > 0 on this interval where 

1-Z(0) + 0Z^(0) 
(1-Z(0))^ 

The denominator is always positive. The numerator is positive if Z^(0) > 0. For 
a standard Type 1 PRC (with a single local extremum), this will occur if 0 is large 
enough (i.e., larger than the minimum point of the PRC; see Fig. lb). For our choice 
of parameters this occurs when 0 > 0.75 (Fig. lb) where the PRC is increasing. On 
the remaining interval, the expression 1 — Z(0) is >1. So if Z^(0) > —1/0 > —4/3 
on [0, 0.75], then /^(0) would also be positive and / could then be inverted on [0, 1] 
(Fig. 4b). However, it is not possible to analytically make this estimate since we have 
no closed form expression for Z(0). We confirmed numerically though that Z^(0) > 
—4/3 in this interval, hence f\<p) is positive on [0, 1]. Therefore, the function / 
can be inverted on [0, 1]. The numerically obtained inverse function /~Ms shown 
in Fig. 4b. Hence, the phase of cell A (3.1a) in cycle n -\-\ can be obtained from its 
value in cycle n from 

0,+i = f{n{f-\4>n))) ^ n(4>n) (3.13) 

In general, the function / (3.12) and the map 77 (3.13) can be defined for networks 
consisting of either identical or non-identical neurons. Here we have considered only 
the networks of identical neurons in this section. The generalization to networks of 
non-identical neurons is considered below in Sect. 3.3.3. 

We can now assess the existence and stability of fixed points of the maps (3.10) 
and (3.13). We numerically solved the map (3.10) using MATLAB to predict the 
activity of two identical M-L neurons coupled with reciprocal inhibition. We also 
numerically solved the differential equations governing the activity of the neurons 
using XPPAUT [20] . We let the maximal synaptic conductance g equal 0. 1 and use 
the PRCs of the neurons obtained for this value of synaptic strength. We first find the 
fixed points of the map by solving the fixed point equation (3.11). The two sides of 
Eq. (3.11) are plotted in Fig. 4a. They intersect only at one point 0* = 0.598, which 
corresponds to the intrinsic phase of cell A (3.2a) at the steady state. The firing period 
of cell A can be obtained from Eq. (3.3a) evaluated at this intrinsic phase. This value 
is also equal to the period of B and will be referred to as the period of the coupled 
network (Pst). The activity phase 0* of cell A (3.1a) at the steady state is 0.5 and 
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Fig. 4 Phase locking for static synapses, a The left and right hand sides of the fixed point equation (3.11) 
for two identical neurons. The left hand side (black) is the response of neuron A and the right hand side is 
the response of neuron B at steady state. The intersection gives the fixed point. Note that the black curve is 
the PRC of both neurons, b The relation between the intrinsic phase 0 (3.2a) and the activity phase 
0 (3.1a). c The same graph as panel a plotted as functions of the activity phase 0 using the transformation 
from (p to 4> shown in panel b. d Convergence of the iterates starting with the initial condition 4>o = 0.2 is 
shown in a cobweb diagram. The iterates (in green) converge to the fixed point at the intersection of the 
graph of 4>n+l = ^(4>n) with the line 4>n =4>n+l 



is obtained by using (3.12), corresponding to the anti-phase solution, which agrees 
with the simulations (not shown). In Fig. 4c, the right and left hand sides of the fixed 
point equation (3.11) are plotted as functions of the activity phase using (3.12). They 
intersect at 0* = 0.5. In Fig. 4d, we show the cobweb diagram for the map (3.13), 
starting with the initial condition 0o = 0-2 leading to convergence to the stable steady 
state of 0* = 0.5. For this case, the system locks in the anti-phase state because the 
two neurons and the two synaptic strengths are identical. 

3.2 Maps Using Dynamic Synapses or Steady-State Synaptic Plasticity Profiles in 
One Synapse 

In this section we derive maps to predict the network activity in the presence of 
synaptic plasticity. We now let the synaptic strength from cell A to cell B be constant 
and the strength from cell B to cell A exhibit plasticity. 
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0 200 400 600 800 1000 1200 1400 
time (msec) 



Fig. 5 Two-cell network with synaptic plasticity in one synapse, a Voltage traces obtained from simu- 
lations of the M-L neurons when the A to B synapse is of fixed strength and B to A synapse changes 
according to the plasticity model (2.3). c The evolution of the plasticity variables r, w, r„, and Un accord- 
ing to the activity of neuron B. d Voltage traces obtained from simulations of the M-L neurons when the 
A to B synapse is of fixed strength and B to A synapse changes according to the steady-state plasticity 
profiles given by (2.4). b & e Network connectivity diagram corresponding to the simulations shown in a 
& d. The parameter values for the plasticity variables are rj =2,t2 = 190, 13 = 2, T4 = 190 

The correct method for deriving the map is to assume that the strength of the 
synapse from B to A changes according to plasticity dynamics given in (2.3). How- 
ever, often in experiments it is easy to measure the steady- state response of a synapse 
at different input frequencies without knowing what the underlying dynamics are 
that give rise to this steady-state value. That is, it is possible to measure the steady- 
state synaptic plasticity profile g{P) obtained from Eq. (2.5). We therefore consider 
two different approaches in the derivation of the map. In the first derivation we as- 
sume that the strength of the B to A synapse is determined by the plasticity dynam- 
ics given in (2.3), whereas, in the second approach, we assume that the strength of 
this synapse obeys the steady-state synaptic plasticity profile g-^iP) given by (2.5) 
(Figs. 5b and 5e). The first approach allows the transients due to different initial con- 
ditions to potentially play a role in the convergence of the map to a fixed point. We 
show, however, that both approaches produce the same result. 

When plasticity is included in the B to A synapse, the synaptic strength is no longer 
constant. Hence we cannot use a unique PRC for neuron A. Instead, we define a PRC 
as a function of two variables, where the phase at which the synapse is received and 
the strength of the synapse determine the response of the neuron. We denote this by 
Za(0, g). The PRC of neuron B is obtained for a constant synaptic strength ^a^b 
and is denoted by Z^{6). 
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We will now determine the phase of neuron A and the network period for the two 
models where the B to A synapse either 

i. changes according to the dynamics of the plasticity variables r and u and is given 

by gB^ArnUn.OY, 

ii. obeys the steady-state synaptic plasticity profile gB(P) = ^B^A^max(^)^min(^). 

We start with the derivation of the map using the dynamics of plasticity variables 
(case i). The voltage traces of the neurons A and B and the evolution of the plasticity 
variables of neuron B obtained from simulations are shown in Figs. 5a and 5c, re- 
spectively. In this case, the response of neuron A in cycle n depends on the values of 
the plasticity variables in this cycle. Assume that we know the values 0^, and Un. 
Then we can compute the period of neuron A in cycle n using the expression 



(3.14) 



We can next modify Eq. (3.5) by rewriting as given in (3.14) to obtain the 
intrinsic phase of neuron B (3.2b) in cycle n as 



Go 



(l - Zp,{(t)n,gB^PJnUn) (pn) 



The equation giving the cycle length of neuron B becomes 



Qn = Go(^l - ^b(^|^(1 - ZA(0n,^B^Arni/n) -0n))) (3.15) 

in cycle n. Using Eq. (3.7) together with the above equations gives a 3D map for the 
evolution of the intrinsic phase of cell A (3.2a) and the synaptic plasticity variables 
from cell B to cell A 



Pn+l 



rn+\ = 



1 - Zb^-^(1 - ZA(0n,^B^A^nW^) -0^)^ 
- 1 + ZA(0n, gB^PJnUn) + 0^ 

\-{\-rne-'^''') 

^^(^^^ ^A(0^,<gB^Ar^i^n) -0n)^^ 



X exp 



-tn 



^2 



(3.16) 



Un+l = U -{U -l^{l-Un)e-'^''') 



X exp 



Go(^l - ^^{^^^ ~ ^A(0^,<gB^A^/|i^n) -0n)^^ 
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The first equation is the same as (3.7) except that now Za is a function of two 
arguments. The second and third equations are computed using (2.3) over one cycle. 
The compHcated expressions in the exponential of both equations are the time Qn — ta 
recast in terms of (l)n,rn, Un where Qn is given in Eq. (3.15). 

We next derive the map for case ii where the synapse from neuron A to neuron B 
has a constant strength at each cycle, while the synaptic strength from neuron B to 
A changes according to the steady-state plasticity function ^b(-^)- The voltage traces 
of the neurons A and B obtained from simulations are shown in Fig. 5d. In this case, 
instead of the depression and facilitation variables, we can use the cycle length of 
one of the neurons to derive the activity map. We assume that we know the values 0^ 
and Pn- Then the intrinsic phase of neuron B (3.2b) in cycle n can be found by using 
(3.4a) as 

On = {Pn-(l>nPQ)IQQ (3.17) 

Plugging this into (3.3b) immediately yields the expression for the cycle length of 
neuron B in cycle n as 

Qn = Qo[l - ZsiiPn - (pnP0)/Q0)] (3.18) 

We can now obtain the intrinsic phase of neuron A (3.2a) in cycle n + 1 using 
Eq. (3.4b) as 

0;.+l = (Qn - dTn)/Po = (Qn " On Qo) / Pq (3.19) 

We use this phase to obtain the cycle length of neuron A in cycle + 1 as 

Pn+l = Po[l - ZA(0n+l,^B(Gn))] (3.20) 

Similar to Eq. (3.14), the period of neuron A is determined by Za which is a 
function of two variables. However, in this case the synaptic strength received by 
neuron A in cycle n -\- 1 depends directly on the cycle length of neuron B in cycle n. 

The map for the activity of the network can be obtained by plugging (3.17) and 
(3.18) into (3. 19) and (3.20) as 



-5^+0n, (3.21) 



Hence, the map (3.16) is reduced to a 2D map for the intrinsic phase and cycle 
length of neuron A. A fixed point (0*, r*, w*) of the 3D map (3.16) corresponds to a 
1 : 1 solution. This 1 : 1 solution is also represented by a fixed point of the 2D map 
(3.21) which occurs at (0*, P*), where is the steady-state value obtained from 
(3.17)at (0*,r*,w*). 

To assess numerically the existence and stability of the fixed points of both the 
2D map (3.21) and the 3D map (3.16), consider two identical neurons coupled with 
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Fig. 6 A comparison of the ID (3.7), 2D (3.21) and 3D (3.16) maps, a The steady-state phase of the 
neuron A, (/)st from map (3.7), ^jyn from map (3.16), (/)ss from map (3.21), shown as a function of the 
intrinsic period of both neurons (changed simultaneously), b The network period as a function of intrinsic 
periods corresponding to the same maps, c The relation between the network period and the phase of A for 
the same maps. The phase of A reaches a minimum at the network period equal to the preferred period of 
neuron B. The results of the two maps with plasticity ((3.14) and (3.19)) overlap in all panels 



asymmetric synapses. Let the synaptic strength from neuron A to B be fixed at 
^A^B =0.1. We use parameters for the plasticity variables that yield the steady- state 
plasticity function ^b(^) with a peak at the period 169.5, as shown in Fig. 2. Denote 
the steady- state network period and phase of neuron A from the 3D map (case i) as 
Pdyn and 0dyn, respectively, and the corresponding values from the 2D map (case ii) 
as Pss and 0ss- Similarly, for static coupling, denote the steady-state network period 
as Pst and phase of neuron A as 0st- 

Figure 6 shows the steady- state phase of neuron A and the network period obtained 
from the ID map (3.7), the 3D map (3.16) and the 2D map (3.21), for a set of intrinsic 
periods (varied simultaneously in both cells). In Fig. 6a, the steady- state phase of 
neuron A is plotted as a function of Pq- The maps with plasticity (cases i and ii) yield 
the same steady-state phase of neuron A; this phase is not constant but is a function of 
the intrinsic period (green and black), in contrast to the static case where the network 
always has an anti-phase solution (dashed red line). This variation in phase depends 
on the values of the steady-state plasticity profile g^{P) (further explained below). 
Figure 6b compares the steady-state network period obtained from the three maps. 
The periods obtained from the maps with plasticity are again the same and they are 
slightly different from the periods obtained from the static map. The blue dashed line 
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is the Po = ^network line. The network period is always larger than the intrinsic period 
in all cases, due to the selection of the PRC (that the inhibitory input always delays 
the next firing time). Figure 6c relates the steady-state phase of neuron A with the 
network period. 

We now examine how the steady- state phase of neuron A changes with respect 
to changes in the intrinsic period. The phase of neuron A depends on the value of 
the synaptic strength received from neuron B at the steady state. This value is deter- 
mined by 2*, the steady-state firing period of neuron B, which equals the steady-state 
network period P*. When this value equals ^a^b =0.1, then anti-phase solutions 
occur. This happens for two sets of coupled neurons, where the red dashed line inter- 
sects green and black curves (Figs. 6a and 6c). Between these two points, the synaptic 
strength received by neuron A, given by (G*), is larger than ^a^b • Since the cells 
are identical, the neurons must give equal amount of response (so that their steady- 
state firing periods will be equal) for a steady-state solution to occur. When both 
synaptic strengths are equal, both neurons have steady-state phase at 0.5. However, if 
^B (G*) > ^A^B, then neuron A receives stronger synaptic input than neuron B. This 
difference can be balanced if neuron A receives this synaptic input at a phase that 
yields less response. As the PRCs of the neurons are decreasing around the phase 
0.5, neuron A needs to phase lock at a phase smaller than 0.5. This explains why 
phase of neuron A decreases between these intersection points. A similar argument 
holds when ^b(G*) < ^a^b- 

The phase of neuron A reaches a minimum when the synaptic strength reaches 
a maximum. As can be seen in Fig. 2, the synaptic plasticity profile has its peak at 
169.5. Therefore, the minimum phase of neuron A is observed at the network period 
169.5 (Fig. 6c). The network period of 169.5 is obtained when two cells with intrinsic 
periods 141.8 are coupled (Fig. 6b). 

3.3 Maps Using Steady-State Synaptic Plasticity Profiles in Both Directions 

Let both reciprocal synapses have short-term plasticity. The map involving the synap- 
tic plasticity variables (2.3) that generalizes (3.16) would now be 5D. But given the 
results from the previous section showing that the simplified map using the steady- 
state synaptic plasticity profiles provides the same stable output, we derive only the 
2D map associated with the latter. We again start with the intrinsic phase 0^ (3.2a) 
and cycle length of neuron A in cycle n. Equation (3.17) can still be used to obtain 
the intrinsic phase of neuron B (3.2b), On, in cycle n. However, the cycle length of 
neuron B is now given by the equation 



in cycle n, since the synapse from neuron A to B also has plasticity and depends on 
Pn . The cycle length P and intrinsic phase <p of neuron A in cycle + 1 is given by 



Qn = Qo[l-Z^{en,gK{Pn))] 



(3.22) 



1 
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Pn+l=n2{(t)n, Pn) = Po[l " ZA(0n+l , (G^))] 



(3.23) 



■Po 
gB( Go 



1-Zb( —(Pn- P0(t>n),gA(Pn) 



1-Zb( —(Pn-Po4>n),gA(Pn) 



Equation (3.23) determines the values of P and 0 when both synapses have plas- 
ticity. In the case where the two cells are identical, Za(-) = Zb(-) = Z, this map 
simplifies to 

1 

0^+1 = /7i(0^, Pn) = —(Qn -Pn + PqM 

Po 

= l-z(^^(Pn- PoM, gA(Pn)^ " ^ + 
Pn+1 = n2{(t)n. Pn) = Po[l " Z(0^+i , (G^))] 

= Po 



(3.24) 



1-Z( 1-Z( -l(P^-Po0n),^A(/^^)j-^+( 



gBiPo 



1-Z( ^-0„^A(/^n) 



We now explore whether these equations yield stable fixed points and, if so, how 
changes in the synaptic profiles affect the resulting phase- and period-locking of the 
network. 

To be able to have explicit control of the preferred frequency of the synapses, 
instead of using Eq. (2.5) for g(P), we assume that the steady-state synaptic pro- 
files obey Gaussian functions gA(') (for the A to B synapse) and ^b(') (for the B 
to A synapse) (2.6) with peaks (preferred periods) Pa and Pb, respectively. Equa- 
tions (3.24) define two surfaces TJiicpn, Pn) and 772(0^, Pn) which can be plotted in 

. We plot two 3D coordinate systems to be able to visualize the evolution of the 
2D map. We show three iterations of the map (3.24) in Fig. 7. The values (cpn, Pn) 
in cycle n are located on the x-y axes. These values are mapped through the sur- 
faces Pn-\-i = n2((pn, Pn) (Fig- 7a) and 0^+i = TJiicpn, Pn) (Fig- 7b) to the next it- 
eration points (0^+1, Pn-\-i) in cycle n -\- I. Start with the initial condition (0o, Po) 
which is shown in both coordinate systems. The image of (0o, Po) on the surface 
0^+1 = 77i (0^ , Pn) gives the next intrinsic phase value 0i , and the image of (0o, Po) 
on the surface P^+i = 772 (0^, Pn) gives the next cycle length Pi (shown by the ver- 
tical lines with one arrow). These 0i and Pi values are located, respectively, on the 
X and 3; axes of both coordinate systems (shown by the inclined lines with one ar- 
row). The point (0i, Pi) is then located on the x-y axes in both coordinate sys- 
tems and mapped to the point (02, P2) by the same procedure (shown by the lines 
with two arrows). We are able to geometrically observe the iterations (only three 
shown) approach a fixed point; hence this is a generalization of cobwebbing for the 
2D map. 
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Fig. 7 Cob webbing diagram of the 2D map (3.24) for two identical cells (Pq = Qq) and distinct synaptic 
plasticity profiles (Pa = 150, Pq = 190) shown in two coordinate systems. The period Pi and the intrinsic 
phase 01 of neuron A in cycle 1 is obtained by evaluating the initial condition {(f)Q, Pq) on the period 
surface P^+i = 772 (0n? Pn) (a) and the phase surface (pn+l = (4>n, Pn) (b)- The point (01, Pi) is then 
projected back to the x-y axis in both coordinate systems and mapped to the point (02 ^ Pi) with the same 
procedure. Lines with one arrow correspond to the first and lines with two arrows correspond to the second 
iteration 



The fixed point equations of the map (3.24) in a 1 : 1 firing condition are 



1-Zb 
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Go 



Go 
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Go 



Go 



Po 
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(3.25) 



These simpHfy to 



Po 



(3.26) 



P* = Po[1-Z(0*,^b(P*))] 



for identical cells. 

The fixed point of this 2D map occurs when (pn = (pn-\-l and Pn = Pn-\-i- We can 
visualize how the fixed point is obtained. For this purpose, we plot the surfaces for the 
evolution of intrinsic phase and period (previously drawn on separate coordinate axes 
in Fig. 7) on the same coordinate axis, above and below the z = 0 plane, and denote by 
the axes zi and Z2, respectively in Fig. 8. The equality 0^ = 0^+i is satisfied when the 
surface zi = ni(x,y) and the plane z\—x intersect. Denote this intersection curve 
as Ci. Similarly, the equality Fn — Pn-\-i is satisfied when the surface Z2 = n2(x, y) 
intersects the plane Z2 = y (denoted as C2). The curves Ci and C2 are shown in black 
above and below the z = 0 plane. The fixed point of the map lies on both curves; 
hence it lies on the intersection of Ci and C2. The projections of Ci and C2 on the 
z = 0 plane are shown in the figure together with the iterations (red dots) approaching 
the fixed point at their intersection. 
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Fig. 8 Fixed points of 2D (3.24) map when Pq = Qq obtained by solving (3.26). The surfaces for the 
evolution of period and intrinsic phase of the 2D map with synaptic preferred periods Pa = 150, Pb = 190 
are drawn above and below the z = 0 plane denoted by the axes zi = Pn+\ and Z2 = 4>n+\^ respectively. 
The equality Pn = Pn+l is satisfied when the surface zi = n2{x, y) (colored surface on top) and the 
plane zi = y (gray-scaled plane on top) intersect. Similarly, the equality (pn = cp^^i is satisfied when 
the surface Z2 = ni(x, y) (colored surface on bottom) intersects the plane Z2 = x (gray-scaled plane on 
bottom). These intersections yield the two black curves above and below the z = 0 plane. The fixed point 
of the map lies on the intersection of the two fixed point curves. The projections of these curves on the 
z = 0 plane are shown together with the iterates (red dots) approaching the fixed point at their intersection 
in the order enumerated in the figure 

The stability of the fixed point can be examined using the Jacobian of the 2D map 
(3.24). If the eigenvalues of the Jacobian at the fixed point are located inside the unit 
circle, the fixed point is stable. For our choice of parameter values, the fixed point 
can be shown to be stable. 

3.3. 1 Phase and Period Locking for Different Synaptic Plasticity Profiles 

Having determined a method for calculating the steady- state network period and 
phase, we now determine how these quantities depend on various network param- 
eters. For simplicity, in this section we consider identical neurons. We use the 2D 
map (3.24) to obtain the network phase and period when both synapses have plas- 
ticity. For comparison, we also obtain the same from the ID map (3.10), when the 
synaptic strength is fixed. 

We are interested in how differences in the plasticity profiles of the two synapses 
affect the network period and phase of neuron A (Figs. 9a 1 and 9b 1). The distinct 
plasticity profiles (Fig. 9al) are produced by simply shifting one profile along the 
intrinsic period axis. In the non-identical case, the plasticity profiles are chosen to 
approach the same value at the tails (Fig. 9al) and, therefore, for small (and large) in- 
trinsic periods, 0ss = 0-5 due to identical synaptic strengths (Fig. 9a3). As the intrin- 
sic period is increased, the difference between gA(P) and ^b(^) first increases until 
P = Pjs^ and then decreases to zero when P = Pgq (Fig- 9al). This causes 0ss to in- 
crease from 0.5 to 0.58 until Pss = and then decrease to 0.5 again when Pss = ^eq 
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Fig. 9 Period and phase locking when both synapses follow the synaptic plasticity profile. Dashed line 
in all panels shows the case with two static synapses, al Synaptic plasticity profiles of the two synapses 
chosen to have different preferred periods at 150 and 190. a2 Network period as a function of the intrinsic 
periods. a3 Phase 0 of neuron A with respect to B as a function of intrinsic period. bl-b3 Same as al-a3 
but with identical synaptic plasticity profiles (preferred period at 170) 



(Fig. 9a3), since the weaker synapse from B to A is balanced by a phase that yields 
more response (more detail is explained in Sect. 3.2). For firing periods greater than 
Peq, the opposite relation holds, causing 0ss first to decrease to 0.41 and then increase 
back to 0.5. In contrast to 0ss varying between 0.41 and 0.58, 0st is always fixed at 0.5 
due to identical neurons and synapses. Since the values of the plasticity profiles at the 
tails are less than the strength g = 0. 1 of the static synapses, Pss is slightly smaller 
than Pst for small (and large) intrinsic periods (Fig. 9a2). For a range of intermediate 
intrinsic periods, when the network synapses have plasticity, Pss is almost equal to 
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the network period with static synapses Pst (Fig. 9a2). The balancing effects of the 
two synaptic profiles (gA(P) being greater, ^b(^) being smaller than g for Pss < ^eq 
and gA(P) being smaller, gB(P) being greater than g for Pss < ^eq) causes Pss and 
Pst to be almost equal for intermediate intrinsic periods. Thus, this choice of synaptic 
plasticity profiles provides the network the ability to produce a range of distinct phase 
relationships as the intrinsic period changes (Fig. 9a3). Note that the steady-state net- 
work period remains almost equal to its value as if no plasticity is included (Fig. 9a2). 

In the case of identical plasticity profiles, the neurons have the same preferred 
periods and the values of the plasticity profiles again approach 0.075 at the tails 
(Fig. 9b 1). This causes Pss to be smaller than Pst for small and large intrinsic pe- 
riods (Fig. 9b2). For intermediate firing periods, the opposite holds. In contrast to 
the almost linear change in Pst, Pss changes nonlinearly as a function of the intrin- 
sic periods. Also, in contrast to the nonlinear change in Pss, the phase of neuron 
A is fixed at 0.5, because both the neurons and their plasticity profiles are identical 
(Fig. 9b3). Hence, depending on the choice of plasticity profiles, the network cou- 
pled with synaptic plasticity can have the same period but different relative phases 
(Fig. 9al-a3), or the same phases but different periods compared to the network cou- 
pled with static synapses (Fig. 9bl-b3). 

3.3.2 Conditions for Phase or Period Constancy 

Short-term synaptic plasticity profiles are subject to change by neuromodulation and 
other long-term modifications [21]. In the previous section, we showed that as the 
synaptic plasticity profile changes, the network can maintain the network period or 
the relative activity phases among the network neurons. In this section, we examine 
the conditions on the steady-state synaptic plasticity profiles that would allow the 
network to maintain either a constant period or a constant phase. 

For this purpose, we make use of the fixed point equations for identical cells (3.26) 
obtained from the 2D map. The phase 0* in (3.26) stand for the intrinsic phase of neu- 
ron A (3.2a). We use Eq. (3.12) and rewrite (3.26) as implicit functions of the steady- 
state phase of A 0, network period P and synaptic preferred periods Pa and Pb as 



Fi(Pa,Pb,(I>,P) = P-Po 
F2(Pa,Pb,^,P) = P-Po 



, P -(l)P 

1-Z( ,gA(P) 



1-Zl ^^gB(P) 



Let P(Pa, Pb, 0, P) = (Fi(Pa, ^b, 0, P), FiiPA, ^b, 0, ^)). At the fixed point, 
F(Pl, P*, 0*, P*) = (0, 0). We would like to solve this equation for Pa and Pb as a 
function of P and 0. Using the Implicit Function Theorem, the condition that needs 
to be satisfied is dQt(Dp^^p^F) 7^ 0 at (P^, Pg , 0*, ^*) where 



■ dFi dFi ■ 

dPA dPB 

dF2 dF2 

.dPA SPb. 



(3.27) 



,P*) 



The function Fi does not depend on Pb, hence 9Pi/9Pb =0. So, for the deter- 
minant to be nonzero, both 9 Pi /9 Pa and 9P2/9Pb have to be nonzero. These terms 
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are given as 



3^1 
dF2 



= ^0^ n ^S^V ) 



One condition for the determinant to be nonzero is dZ/dy(x, y)\ .p* 7^ 0; 

^ A' B ' ^ 

that is, the response of the neuron to perturbations should change with the change 
in the strength of the perturbation. This is a standard assumption on phase re- 
sponse curves with small perturbations. The other two conditions to be satisfied 
are 9^a/9^aI(p*,p*,0*,p*) 7^ 0 and 9<gB/9^B 7^ 0, which, upon using 
Eq. (2.6), are equivalent to Pa 7^ and ^ P"^, respectively. In other words, the 
network period should be different from the synaptic preferred periods. 

Under these three conditions, the Implicit Function Theorem guarantees that Pa 
and Pb can be expressed in terms of 0 and P near (P^, Pg, 0*, P*). More pre- 
cisely, there are neighborhoods U of (0*, P*) and W of (P^, Pg) such that, for 
each (0, P) G U, there exists a unique (Pa, Pb) ^ W such that P(Pa, Pb, 0, P) = 
P(Pa(0, P), Pb(0, P), 0, P) = 0. Hence, there is a unique function h = (hi,h2) : 

^ W such that P(/zi(0, P), /z2(0, ^), ^,P) = 0 for every (0, P) g U. 

We can interpret this result in two ways. First, around the fixed point (P^, Pg, 
0*, P*), we can choose (0^ P*) such that P* is fixed and 0^ 7^ 0*, for which there 
exist (PA^ PbO that satisfy the fixed point equations (3.26). Hence, for a specific P*, 
around a point with a phase 0^ there exist synaptic preferred periods Pa' and Pb' 
that enable the network to stay on the level set of P*, while setting the phase equal 
to a new value 0^ In other words, it is possible to keep the network period constant 
and set the network phase to a new value by changing the synaptic plasticity profiles 
of the network neurons. 

The second interpretation is that, around the fixed point (P^, Pg , 0*, P*), we can 
choose a (0*, PO such that 0* is fixed and P^ 7^ P*, and can find (Pas PbO that 
satisfy the fixed point equations (3.26). This enables the network to stay on the level 
set for a specific 0*, while changing the network period to a new value P^ 

In the example demonstrated in Fig. 10, the intrinsic periods of the two neurons 
are kept constant but the two synaptic plasticity profiles are allowed to vary. As be- 
fore, the synaptic plasticity profiles are changed only by shifting them along the pe- 
riod axis. We keep track of different synaptic plasticity profiles by the values of the 
synaptic preferred periods Pa and Pb (the peak of the profile). Figure 10 shows the 
changes in the network period and phase as the synaptic plasticity profiles of the neu- 
rons are varied. The neurons are identical with an intrinsic period Pq of 137. The 
colored curves are subsets of the level sets of the phase; the phase of the network is 
fixed on a curve with a specific color. The gray bands correspond to the level sets of 
the network period. These level sets provide conditions for the network to maintain a 
specific period but have different phase relations, or vice versa, through varying the 
combination of synaptic preferred periods. 
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Preferred Period (^) of A to B synapse 

Fig. 10 Period and phase locking for different steady-state synaptic plasticity profiles. The steady-state 
network period (gray) and phase (colored) are shown as a function of different steady-state synaptic plas- 
ticity profiles. Colored curves correspond to level sets of the phase. The edges of the gray bands correspond 
to the level sets of the network period. The plasticity profile of each synapse is marked by its preferred 
period 

3.3.3 Networks of Non-identical Neurons 

We now examine a network of two non-identical M-L neurons. The neurons are 
chosen to have different intrinsic periods by applying different levels of external cur- 
rent but otherwise using the same parameters. We consider the two cases where the 
synapses are static or they follow steady-state synaptic plasticity profiles and compare 
the predictions of the ID map (3.7) and the 2D map (3.24) with the simulations of the 
corresponding model equations. We let the preferred period of the A to B synapse be 
Pa = 150 and from neuron B to A be Pb = 190 for the case with synaptic plasticity. 
The results are shown in Fig. 11. 

Note that the maps continue to give good predictions when the neurons are not 
necessarily identical. The difference between the simulations (filled circles) and the 
map predictions (open circles) is indistinguishable in most cases. The diagonal corre- 
sponds to coupling of identical neurons. Moving away from the diagonal, the differ- 
ence between the intrinsic periods of the neurons increases and eventually prevents 
the neurons to phase lock in a 1 : 1 manner because the fixed point equation (3.8) is 
not satisfied anymore. These are the limits of the region shown in Fig. 11. Observe 
that the limits determined by the map and the simulations overlap except at one single 
case shown only by an open circle in Figs. 11c and lid. Here, the map predicts that 
a 1 : 1 solution exists, while the simulation does not converge to that. In this case, 
the simulation shows that the firing order between the neurons is not preserved which 
violates the 1 : 1 firing assumption of the map. 
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Fig. 11 Coupling of non-identical M-L neurons. The phase of neuron A (a and c) and the period of the 
network (b and d) for coupled neurons with different intrinsic periods are shown for static synapses (a and 
b; ^ = 0.1) and when the network follows the synaptic plasticity profile (c and d; Pa = 150, Pb = 190). 
The axes are the intrinsic periods of the two neurons. Plasticity adds nonlinearity to the period and phase 
distribution. Filled circles denote simulation results whereas open circles denote the map predictions. The 
map yields predictions very close to the simulations in most cases 



The phase of neuron A equals 0.5 on the diagonal in the static coupling case 
(Fig. 11a). It decreases (resp. increases) linearly as Qo moves down (resp. up) from 
the diagonal. This behavior can be predicted by analyzing Eq. (3.8). In the identi- 
cal network, where = Qq, the activity phases (0* = = 0.5), and the intrinsic 
phases (0* = ^* = 0.598) of the two neurons are equal and hence Za(0*) = Zb(^*). 
If the solution is perturbed such that Pq > Qo, then the response of neuron A to 
synaptic inputs from neuron B must be smaller than the response of neuron B for the 
Eq. (3.8) to be satisfied. The PRC of the neurons has a negative slope at this intrinsic 
phase 0* (Fig. lb). So, the intrinsic phase 0 of neuron A in the perturbed solution 
must be smaller than 0* for Za(0) to be smaller than Za(0*). As the function (3.12) 
relating 0 and 0 is monotone increasing, the activity phase 0 of neuron A in the 
perturbed solution must also be smaller than 0*. Hence, as the difference Pq — Qo 
increases (resp. decreases), the phase of neuron A decreases (resp. increases). The 
period of the network increases linearly as the intrinsic periods increase in the static 
coupling case (Fig. lib). Due to symmetry in the synaptic strengths, the distribution 
of the period is symmetric with respect to the diagonal. 
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When the synapses are plastic, some 1 : 1 phase-locked solutions that existed with 
static coupling no longer exist, while new solutions may emerge (Figs. 11c and lid). 
Due to asymmetry in the synaptic plasticity profiles, the upper bound for the dif- 
ference in intrinsic periods that allow a 1 : 1 phase-locked solution varies. This can 
be seen by comparing the circles in the top row and rightmost column of Figs. 11c 
and lid. At the right top corner, Pq = Go = 181, and the network has an anti-phase 
solution. If 2o is fixed, while Pq decreases, the network continues to phase lock in a 
1 : 1 solution for Pq > 152.1. On the other hand, if Pq is fixed, while Qq decreases, 
then the network phase locks in a 1 : 1 solution only when Qq > 174.8. Although 
the absolute difference between the intrinsic periods are equal, different plasticity 
profiles causes convergence in one case but not the other. This can be understood 
by considering (3.25). For the identical cell case where = Go = 181, the net- 
work period is equal to P* = 219.5. Due to the selection of the plasticity profiles, 
^a(^*) < ^b(^*), since P* is close to Pr = 190 than it is to Pa = 150. As a re- 
sult, neuron A receives stronger synaptic input from neuron B at the steady state (as 
^b(^*) determines ^b^a). The firing periods of both neurons must be equal at the 
fixed point. This is only possible if neuron B receives synaptic input at a phase that 
yields a larger response than that of neuron A. Hence, although the neurons are identi- 
cal, the difference in their plasticity profiles causes a phase-locking solution different 
from anti-phase. Assume now that Go > ^0- Then the relation gA(P) < <^b(^) will 
still hold as P will stay close to P*. In this case, the synaptic strength received by 
neuron A will be larger, while its intrinsic period will be smaller than that of B. These 
two opposing effects will let the network continue having a solution until the differ- 
ence between the intrinsic periods are too large to be compensated for and (3.25) 
are not satisfied. On the other hand, if the symmetric solution is perturbed such that 
^0 > Go, then the synaptic strength received by neuron A and its intrinsic period will 
both be larger than those of neuron B. The phase of neuron B must increase further 
and yield a larger response to compensate for these adding effects. But when the PRC 
reaches a maximum in absolute value and starts to decrease, there would be no phase 
value that would compensate for these effects and the network will not be able to have 
a 1 : 1 solution. This explains why the limits of the regions in the case with synaptic 
plasticity are not symmetrical. 

In general, whether (3.25) are satisfied or not depends on the intrinsic periods Pq, 
Go and the values of the PRCs as in the static map case. But in this case the values of 
the PRCs are also determined by two factors, the phase of inhibition received, and its 
strength — which is determined by the network period. Hence, the phase of neuron A 
is a determined both by the interaction of intrinsic periods and the plasticity profiles. 
This is also responsible for the nonlinearity in the distribution of phase. The level 
curves of phase are nonlinear in the case with synaptic plasticity as opposed to the 
linear level curves in the static coupling case. 



4 Discussion 

In the analysis of an oscillatory network, the steady-state activity of the network can 
often be reduced to the study of a return map. The advantage of using maps is that 
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it often allows the network dynamics to be understood by tracking empirically ob- 
servable characteristics such as period and phase. Here, we derive such a map for 
a two-cell network coupled with inhibitory synapses with the goal of understanding 
how short-term synaptic plasticity and other factors determine the network period and 
the relative activity phase of the two neurons. Our results show that the information 
on the network period and phase can be obtained using maps that keep track of ob- 
servable network variables such as the intrinsic periods of the neurons involved, their 
phase response curves and the synaptic plasticity profiles: relationships describing 
how the synaptic strength depends on input frequency. These variables can be read- 
ily determined experimentally with "feed-forward" measurements where the input is 
controlled by the experimenter and the output is measured. For example, the strength 
of a synapse can be measured at all frequencies simply by driving the presynaptic 
neuron at different rates and measuring the postsynaptic current. In fact, the current 
study was motivated by our experimental measurements of these types of network 
variables in the crab stomatogastric pyloric network [22-24]. 

There are several prior works that utilize PRCs and map-based techniques to un- 
derstand phase locking [1-13]. Of particular interest is the result of Cui et al. [5] who 
use a functional PRC (fPRC) that is calculated from actual experimental measure- 
ments of Aplysia pacemaker neurons. Cui et al. show that the fPRC differs from the 
single pulse PRC (as was used in this paper) due to accommodation of the pacemaker 
neurons. They then go on to use the fPRC to study phase-locking in a coupled network 
by deriving a map that encodes how a neuron responds to a period input that arrives 
a fixed time after the firing of the cell. By linearizing about a fixed point of their ID 
map, they find conditions for the existence and stability of 1 : 1 phase-locked solu- 
tions. Their predictions from the fPRC method are better matched to simulations than 
predictions from a conventional single-pulse PRC. Importantly, their fPRC methods 
do not depend on the exact shape of the PRC but rather on the effect on the cycle 
period based on the time the input was given. This is a statistic that is easily found in 
experiments. Moreover, their results are obtained from combining feed-forward pro- 
cesses as opposed to directly studying a feedback map, which they call open-looped 
versus closed-looped. 

Our results complement those of Cui et al. in the sense that we relate cycle-to- 
cycle changes in the period independent of how those changes arise, allowing us to 
also use experimentally obtainable information to derive the maps. Our maps are also 
based on assumptions that are consistent with Cui et al.'s assumption that the closed- 
loop behavior of a system can be predicted by knowing the open-looped behavior of 
some of its components. Our results extend those of Cui et al. and other prior works 
in that we allow the timing of inputs to vary on a cycle by cycle basis that is deter- 
mined by the synaptic plasticity profile of the presynaptic neuron. This results in a 
higher-dimensional map arising by specifically considering the dynamics of synap- 
tic facilitation and depression on a cycle by cycle basis. This yields a 3D map when 
plasticity is present only in one direction of the two-cell network, or a 5D map if 
present in both directions. When we used the steady-state synaptic plasticity profile, 
both cases reduce to a 2D map. For this 2D map, we derived a geometric method that 
generalizes cob webbing in a ID map to allow us to study the existence and stability 
of fixed points. For a generic ID map, 77 (x), the intersection of the curve y = n(x) 
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with the curve y = x, and the slope at that point, determine existence and stabiHty of 
the fixed point. In our generaUzed 2D case, given maps /7i(x, y) and 772 (x, y), it is 
the intersection of these surfaces with appropriate planes that yield two curves. It is 
the intersection of the projection of these two curves onto a common plane that deter- 
mines existence of the fixed point. Stability is more complicated than just checking 
the slopes at the point of intersection. We showed how it could depend on both the 
PRC and the synaptic plasticity profile. 

In this study, we considered a general form of short-term synaptic plasticity 
which is a combination of short-term facilitation and depression. We modeled such 
a synapse using an ad hoc model as described previously [16]. The advantage of this 
model is that the extent to which facilitation or depression is a dominant factor can 
be simply determined by changing the model parameters. Our analysis progressed 
through a network of two neurons with static synapses, the same network but with 
one synapse having plasticity and finally with both synapses showing plasticity. The 
analysis of a two-cell network with static synapses yields a ID map [6, 8]. Including 
synaptic plasticity increases the dimension of the map because variables underlying 
synaptic dynamics must be tracked as well. The change in synaptic strength due to 
the plasticity means that the PRCs of the neurons also change. Our analysis shows 
that these higher-dimensional maps can accurately predict the steady- state phase and 
period of the network, as seen in comparisons with numerical simulations of the un- 
derlying ODEs. 

In experimental measurements, synaptic plasticity profiles are often measured us- 
ing repetitive input pulses or waveforms and reported at steady state, i.e., the steady- 
state strength of the synapse is known for each stimulation frequency [23, 25, 26]. 
In most cases, the mechanisms that underlie these synaptic dynamics are unknown 
and it is therefore impossible to track how synaptic strength changes as a function of 
frequency on a cycle-to-cycle basis. One of the interesting findings from our work is 
that the prediction of the higher-dimensional map obtained when using dynamics of 
the synapse is the same as a lower-dimensional map that uses only the steady-state 
plasticity profile. In other words, the network output is dependent on the steady- state 
strength independent of the mechanisms through which this synaptic strength is ac- 
tually generated. In turn, this allows an experimentalist to understand the effects of, 
say a synaptic neuromodulator, on the network output simply by understanding the 
effect on a single component such as the synaptic plasticity profile. 

The results of our maps help us understand the role of synaptic dynamics in deter- 
mining the relative phase between two neurons in an oscillatory network. For exam- 
ple, neurons in the crustacean pyloric oscillatory network, involve multiple recipro- 
cally inhibitory connections. Pyloric oscillations are quite stable in individual prepa- 
rations and are generated by a pacemaker group of neurons (AB/PD) which make 
reciprocally inhibitory connections with a single follower neuron, LP. The analysis of 
this reciprocally inhibitory network provided the motivation for the current study. As 
in other rhythmic motor networks, the pyloric network neurons maintain a constant 
phase relationship even when these phases are measured in different animals [27]. 
Surprisingly this tight phase relationship is maintained despite a large variability in 
the pyloric cycle period (1-2.5 s) across preparations. In fact, different preparations 
differ both in the intrinsic periods of the neurons involved as well as the synaptic 
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plasticity profiles. The results of the current study indicate that the pyloric network 
could maintain constant phase relationships, even in different animals, by tuning the 
synaptic plasticity profiles along the level sets of phase (Fig. 10). Alternatively, if the 
relative activity phases of the neurons involved in producing the network oscillations 
are not an essential component of the network output, but the network must maintain 
a constant period, the maps we have derived can be used to establish the relationships 
that could produce a constant frequency output. These are plausible strategies for all 
rhythmic motor networks in which the output is tightly constrained by the proper 
phase of muscle movements to produce meaningful behavior. 

An interesting implication of our results is that if the network period coincides 
with the synaptic preferred periods, it is not possible to uniquely prescribe the synap- 
tic profiles in terms of the network period and the relative phase of the neurons 
(Eq. (3.27)). If the level sets of phase, described in Fig. 10, provide a unique rule 
for the network to tune its synaptic plasticity profiles for phase maintenance, then the 
network period should avoid the synaptic preferred periods. Additionally, by avoid- 
ing the periods at which the synaptic strengths are maximal, the network can operate 
with a larger degree of flexibility and perhaps more efficiently. This is in fact the case 
for the synapses between the AB/PD pacemaker neurons and the follower LP neuron 
in the crustacean pyloric network. The network period is around 1-2.5 s, in a range of 
values that is larger than the preferred periods of the synapses (~0.5 Hz) [23]. Hence, 
our findings give an insight for this experimentally observed fact. 

In conclusion, we have shown that the frequency-dependent information on 
synapses can be combined with the PRCs of oscillatory neurons to predict the ac- 
tivity period and phases of a coupled network using maps derived from empirically 
observable relationships. It is plausible that a similar approach can be used whenever 
there is frequency-dependent information about the network components to construct 
maps that predict the activity of an oscillatory network, even when the synapses in- 
clude excitatory connections or obey different plasticity profiles. In relationship to 
the crustacean pyloric network that motivated this study, current experimental work 
in our lab involves measuring the changes in the synaptic plasticity profiles and the 
neuronal PRCs in the presence of different neuromodulators to see whether the maps 
derived here can predict how the network output changes in the presence of these 
modulators. 
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